\chapter{Polarizable embedding calculations}\label{ch:embedding}

\index{polarizable embedding}\index{electrostatic embedding}\index{PE}\index{environment model}\index{embedding model}\index{multiscale modeling}\index{embedding}\index{QM/MM}\index{quantum mechanics / molecular mechanics}\index{solvent effects}\index{enviroment effects}\index{polarizable density embedding}\index{PDE}
This chapter introduces the polarizable embedding (PE)\cite{pemodel1,pemodel2} and
polarizable density embedding (PDE)\cite{pde1,pde2,pde3} models as implemented
in the PE library~\cite{pelib} included in the \latestrelease\ release.
Methods available are: PE-HF~\cite{pescf}, PE-DFT~\cite{pescf},
PE-MP2/SOPPA~\cite{pesoppa}, PE-MCSCF~\cite{pemcscf}, PE-MC-srDFT~\cite{pesrdft}
and PE-CC~\cite{pecc}.
The PDE model can be used with all of the same methods (see subsection~\ref{subsec:pde})~\cite{pde1,pde2,pde3,pdecc}.
The implementation uses the Gen1Int library to calculate one-electron
integrals~\cite{gen1int} which is also included in the \latestrelease\ release. The
first section gives some general considerations about the model and
implementation. In the second section we introduce the input format using
basic examples.
We also refer to our tutorial review on the use of polarizable embedding for modeling of response
properties of embedded molecules~\cite{petutorial} (see \url{https://doi.org/10.1002/qua.25717}
or alternatively \url{https://arxiv.org/abs/1804.03598}).

\begin{center}
\fbox{
\parbox[h][\height][l]{12cm}{
\small
\noindent
{\bf Reference literature:}
\begin{list}{}{}
\item Polarizable embedding (PE): J.~M.~Olsen, K.~Aidas and J.~Kongsted, \newblock {\em J.~Chem.~Theory~Comput.}, {\bf 6}, 3721 (2010) and J.~M.~H.~Olsen and J.~Kongsted, \newblock {\em Adv. Quantum Chem.}, {\bf 61}, 107 (2011)
\item PE-HF/DFT: J.~M.~Olsen, K.~Aidas and J.~Kongsted, \newblock {\em J.~Chem.~Theory~Comput.}, {\bf 6}, 3721 (2010).
\item PE-MP2/SOPPA: J.~J.~Eriksen, S.~P.~A.~Sauer, K.~V.~Mikkelsen, H.~J.~Aa.~Jensen and J.~Kongsted, \newblock {\em J.~Comp.~Chem.}, {\bf 33}, 2012, (2012).
\item PE-MCSCF: E.~D.~Hedeg\aa{}rd, N.~H.~List, H.~J.~Aa.~Jensen and J.~Kongsted, \newblock {\em J.~Chem.~Phys.}, {\bf 139}, 044101 (2013).
\item PE-CC: K.~Sneskov, T.~Schwabe, J.~Kongsted and O.~Christiansen, \newblock {\em J.~ Chem.~Phys.}, {\bf 134}, 104108, (2011).
\item Damped response (CPP) with PE: M.~N.~Pedersen, E.~D.~Hedeg\aa{}rd, J.~M.~H.~Olsen, J.~Kauczor, P.~Norman and J.~Kongsted, \newblock {\em J.~Chem.~Theory~Comput.}, {\bf 10}, 1164 (2014).
\item Magnetic properties with LAOs: C.~Steinmann, J.~M.~H.~Olsen and J.~Kongsted, \newblock {\em J.~Chem.~Theory~Comput.}, {\bf 10}, 981 (2014).
\item Effective external field (EEF) effects: N.~H.~List, H.~J.~Aa.~Jensen and J.~Kongsted, \newblock{\em Phys.~Chem.~Chem.~Phys.},{\bf 18}, 10070 (2016) and N.~H.~List, PhD thesis, University of Southern Denmark, Odense, Denmark, 2015.
\item Polarizable density embedding (PDE): J.~M.~H.~Olsen, C.~Steinmann, K.~Ruud, and J.~Kongsted, \newblock{\em J.~Phys.~Chem.~A}, {\bf 119}, 5344 (2015), P.~Reinholdt, J.~Kongsted, and J.~M.~H.~Olsen, \newblock {\em J.~Phys.~Chem.~Lett.}, {\bf 8}, 5949 (2017), and P.~Reinholdt, F.~K.~J\o{}rgensen, J.~Kongsted, and J.~M.~H.~Olsen, \newblock {\em J.~Chem.~Theory~Comput.}, {\bf 16}, 5999 (2020)
\item Tutorial review: C.~Steinmann, P.~Reinholdt, M.~S.~N\o{}rby, J.~Kongsted and J.~M.~H.~Olsen, \newblock{\em Int.~J.~Quantum~Chem.}, {\bf 119}, e25717 (2019), \url{https://doi.org/10.1002/qua.25717}.
\end{list}
}}
\end{center}

\section{General considerations}
In {\dalton} it is possible to include the effects from a structured
environment on a core molecular system using the polarizable embedding (PE)
model. The current implementation is a layered QM/MM-type classical embedding model
capable of using advanced potentials that include an electrostatic component as
well as an induction (polarization) component. The effects of the environment
are included through effective operators that contain an embedding potential,
which is a representation of the environment, thereby directly affecting the
molecular properties of the core system. The wave function of the core system
is optimized while taking into account the explicit electrostatic interactions
and induction effects from the environment in a fully self-consistent manner.
The electrostatic and induction components are modeled using Cartesian
multipole moments and anisotropic dipole-dipole polarizabilities, respectively.
The electrostatic part models the permanent charge distribution of the
environment and will polarize the core system, while the induction part also
allows polarization of the environment. The environment response is included in
the effective operator as induced dipoles, which arise due to the electric
fields from the electrons and nuclei in the core system as well as from the
environment itself. It is therefore necessary to recalculate the induced dipoles
according to the changes in the electron density of the core system as they
occur in a wave function optimization. Furthermore, since the induced dipoles
are coupled through the electric fields, it is necessary to solve a set of
coupled linear equations. This can be done using either an iterative or a
direct solver. This also means that we include many-body effects of the total system.

The multipoles and polarizabilities can be obtained in many different ways. It
is possible to use the molecular properties, however, usually
distributed/localized properties are used because of the better convergence of
the multipole expansion. These are typically centered on all atomic sites in
the environment (and sometimes also bond-midpoints), however, the
implementation is general in this sense so they can be placed anywhere in
space. Currently, the PE library supports multipole moments up to fifth order and
anisotropic dipole-dipole polarizabilities are supported. For multipoles up to and
including third order (octopoles) the trace will be removed if present. Note, that
the fourth and fifth order multipole moments are expected to be traceless. In case
polarizabilities are included it might be necessary to use an exclusion list
to ensure that only relevant sites can polarize each other. The format of
the \potinp\ file is demonstrated below.

The PE model is implemented for HF, DFT, MP2, MCSCF, MC-srDFT, and
CC (CC2, CCSD, and CCSDR(3)) wave functions. Singlet linear response may be evaluated
at PE-HF and PE-DFT levels of theory for closed- and open-shell systems,
and at PE-SOPPA, PE-MCSCF, PE-MC-srDFT, PE-CC2, PE-CCSD, and PE-CCSDR(3) levels for
closed-shell systems.
Furthermore, the PE library has been coupled to the complex polarization
propagator (CPP) module to allow damped response properties at
PE-HF and PE-DFT levels of theory for closed- and open-shell systems~\cite{pecpp}.
Triplet linear response properties are
available at the PE-HF and PE-DFT levels for closed-shell systems
only.
Singlet quadratic and cubic response properties are available for PE-HF and PE-DFT
levels of theory for closed-shell systems.
It is furthermore possible to include the effective external field (EEF) effect,
which models the environment polarization induced directly by the presence of an
external field~\cite{peeef,peeef2}.
This leads to properties that are defined in terms of the external field and are
thus comparable to supermolecular calculations.
Magnetic linear response properties using London atomic orbitals (LAOs) are also
available~\cite{pelao}.
In the PE-MP2/SOPPA models, the environment response is taken into account at the
HF/RPA level (see Ref.~\cite{pesoppa} for details).
Note, that point-group symmetry and analytical molecular Hessians are not supported.


\section{Input description}
To include environment effects using the PE model it
is necessary to define the core and environment part of the system: the \molinp\
file specifies the core molecular system and the \potinp\ file, which contains
the embedding potential, defines the environment. Moreover, additional
keywords are needed in the \dalinp\ input file to activate the PE model. To use
default options it is only necessary to include \Key{PELIB} in the \Sec{*DALTON}
section (except for PE-CC which requires additional input (see subsection~\ref{subsec:pecc})).
All other specifications of wave function, properties etc.\ are
unchanged and thus follow the input described in other chapters. For example, to
calculate the PE-HF wave function the following input can be used:
\begin{verbatim}
**DALTON
.RUN WAVE FUNCTIONS
.PELIB
**WAVE FUNCTIONS
.HF
**END OF DALTON
\end{verbatim}
To use non-default options, a \Sec{PELIB} subsection is needed, which should also be placed in the
\Sec{*DALTON} section. For instance, to use the direct solver for induced
dipoles the following input example can be used:
\begin{verbatim}
**DALTON
.RUN WAVE FUNCTIONS
.PELIB
*PELIB
.DIRECT
**WAVE FUNCTIONS
.HF
**END OF DALTON
\end{verbatim}
where the \Key{DIRECT} keyword request the use of a direct solver. See further
input options in Chapter~\ref{ch:general} under the \Sec{PELIB}
input section (subsection~\ref{subsec:pelib}). Furthermore,
Section~\ref{sec:daltoninp} in Chapter~\ref{ch:starting} provides an
introduction to the \dalton\ (and \molinp) input in general. The format of the
\molinp\ file is described in detail in \ref{ch:molinp} and requires no
additional changes to be used in a PE calculation.

\subsection*{The potential input format}
The \potinp\ file is split into three sections: \verb|@COORDINATES|,
\verb|@MULTIPOLES| and \verb|@POLARIZABILITIES|. Additionally, an optional
\verb|@SIMULATION_BOX| section can be included to define the simulation box.
It is required for calculations including periodic boundary conditions.
The format of the \potinp\ file is perhaps best illustrated using an example:
\begin{verbatim}
! two water molecules
@COORDINATES
10
AA
O    -3.328  -0.103  -0.000
H    -2.503   0.413   0.000
H    -4.039   0.546  -0.000
X    -2.916   0.154  -0.000
X    -3.683   0.221  -0.000
O     1.742   2.341  -0.000
H     0.841   1.971  -0.000
H     1.632   3.298   0.004
X     1.291   2.156  -0.000
X     1.687   2.819   0.001
@MULTIPOLES
ORDER 0
6
1    -0.742
2     0.369
3     0.372
6    -0.742
7     0.369
8     0.372
ORDER 1
10
1     0.030    0.328    0.000
2    -0.100   -0.055   -0.000
3     0.091   -0.072    0.000
4    -0.115   -0.109   -0.000
5     0.092   -0.128    0.000
6    -0.284    0.167    0.001
7     0.103    0.049    0.000
8     0.005   -0.116   -0.000
9     0.156    0.028   -0.000
10    0.050   -0.149   -0.000
ORDER 2
10
1    -3.951   -0.056    0.000   -4.577    0.000   -5.020
2    -0.577   -0.053   -0.000   -0.604   -0.000   -0.559
3    -0.558    0.046    0.000   -0.622    0.000   -0.558
4     0.693    0.399    0.000    0.481    0.000    0.233
5     0.549   -0.407    0.000    0.632   -0.000    0.241
6    -4.418    0.280    0.000   -4.112    0.003   -5.020
7    -0.645   -0.004    0.000   -0.536    0.000   -0.559
8    -0.556    0.045    0.000   -0.624   -0.000   -0.558
9     0.930    0.228   -0.000    0.242   -0.000    0.233
10    0.217   -0.166   -0.000    0.964    0.003    0.241
@POLARIZABILITIES
ORDER 1 1
10
1     1.593    0.080   -0.001    2.525    0.001    3.367
2     0.792    0.154    0.000    0.601    0.000    0.592
3     0.720   -0.178    0.000    0.642    0.000    0.575
4     3.497    2.135    0.002    1.845    0.001    1.412
5     2.691   -2.246    0.000    2.554   -0.001    1.429
6     2.282   -0.420   -0.000    1.832   -0.006    3.366
7     0.813    0.138    0.000    0.581   -0.000    0.592
8     0.499   -0.019    0.000    0.861    0.002    0.575
9     4.294    1.269    0.000    1.056   -0.004    1.413
10    0.617   -0.440   -0.000    4.622    0.017    1.430
EXCLISTS
10 5
1  2  3  4  5
2  1  3  4  5
3  1  2  4  5
4  1  2  3  5
5  1  2  3  4
6  7  8  9 10
7  6  8  9 10
8  6  7  9 10
9  6  7  8 10
10 6  7  8  9
\end{verbatim}
Note, that the input is actually not case-sensitive even though we have used
uppercase letters in the example. In the following we describe the
contents of each section.\\

\noindent\texttt{@COORDINATES}\newline
The coordinates section follows the standard XYZ file format so that the
environment can be easily visualized using standard programs. The first line in
gives the total number of sites in the environment and the second line
specifies whether the coordinates are given in \angstrom{} (\verb|AA|) or
\bohr{} (\verb|AU|). The rest of the coordinates section is a list of the sites in
the environment where each line contains the element symbol and x-, y- and
z-coordinates of a site. If a site is not located on an atom, e.g.\ if it is a
bond-midpoint, then the element symbol should be specified as \verb|X|. The
listing also gives an implicit numbering of the sites, so that the first line
is site number one, the second line is site number two and so on. This
numbering is important and used in the following sections.\\

\noindent\texttt{@MULTIPOLES}\newline
The multipoles section is subdivided into the orders of the
multipoles, i.e.\ \verb|ORDER 0| for monopoles/charges, \verb|ORDER 1|
for dipoles and so on.  For each order there is a number specifying
the number of multipoles of that specific order. Note, that this
number does not have to be equal to the total number of sites. This is
followed by a list of multipoles where each line gives the multipole
of a site. The lines begin with a number that specifies which site the
multipole is placed. Only the symmetry-independent Cartesian
multipoles (given in a.u.) should be provided using an ordering such
that the components are stepped from the right, e.g.\
\verb|xx xy xz yy yz zz| or \verb|xxx xxy xxz xyy xyz xzz yyy yyz yzz zzz|.
Note, that the multipoles should in general be traceless, however, for
multipoles up to and including octopoles (i.e.\ \verb|ORDER 3|) the
trace is removed if present. Furthermore, the current implementation is
limited to fifth-order multipoles.\\

\noindent\texttt{@POLARIZABILITIES}\newline
The polarizabilities section is also subdivided into orders,
i.e.\ \verb|ORDER 1 1| for dipole-dipole polarizabilities, which is the only
type supported in the current release. The format is the same as for
multipoles, i.e.\ first line contains number of polarizabilities which is
followed by a list of the polarizabilities using the same ordering as the
multipoles. The polarizabilities should also be given in a.u. In addition,
there is also the exclusion lists (\verb|EXCLISTS|
section). Here the first line gives the number of lists (i.e.\ the number of
lines) and the length of the exclusion lists (i.e.\ the number of entries per
line). The exclusion lists specify the polarization rules. There is a list
attached to each polarizable site that specifies which sites are not allowed
to polarize it, e.g.\ \verb|1 2 3 4 5| means that site number 1 cannot be
polarized by sites 2, 3, 4 and 5.\\

\noindent\texttt{@SIMULATION\_BOX}\newline
The simulation box section is used to define the simulation box type and
parameters. The box can either be orthorhombic:
\begin{verbatim}
@SIMULATION_BOX
ORTHORHOMBIC
a b c
\end{verbatim}
or triclinic:
\begin{verbatim}
@SIMULATION_BOX
TRICLINIC
a b c A B C
\end{verbatim}
Here a, b, and c are the side lengths in \AA{} and A, B, and C are the angles in degrees.

\subsection*{PE-CC example}\label{subsec:pecc}
To run a PE-CC calculation requires an additional \Sec{PECC} section under
\Sec{*WAVE FUNCTIONS} where additional options specfic to PE-CC are placed.
The following input exemplifies the input:
\begin{verbatim}
**DALTON
.RUN WAVEFUNCTION
.PELIB
*PELIB
**WAVE FUNCTIONS
.CC
*CC INP
.CCSD
*PECC
.MXSLIT ! max. no. t/t-bar iterations in solution of coupled t/bar-t eqs.
 200
.MXINIT ! max. no. of steps in the t and t-bar solver, respectively
 4 5
*CCEXCI
.NCCEXCI
 2
*CCLRSD
.DIPOLE
**END OF
\end{verbatim}
For PE-CC calculations
it is obligatory to include the \Sec{PECC} section. Also given in the example
is the maximum number of $t$/$\bar{t}$ iterations in the solution of the
coupled $t$/$\bar{t}$ equations (\Key{MXSLIT}) and the maximum number of
steps in the $t$ and $\bar{t}$ solver (\Key{MXINIT}). For details regarding
the general input for CC calculations we refer to
Chapters~\ref{ch:ccexamples} and ~\ref{ch:CC}.

\subsection*{Pseudopotentials for avoiding electron spill-out}
PE calculations may suffer from electron spill-out, especially when using diffuse
basis functions, due to the lack of non-electrostatic Pauli repulsion from the electrons
of the environment. This can, in some cases, lead to unphysical charge leaks from
the core quantum region to the environment.
A cost-effective solution to prevent this is to use the pseudopotentials developed
in ref. \citenum{pe_pp}.
To add pseudopotentials on classical atoms in the environment, include them in the \molinp file with
\texttt{Basis=pointcharge ECP=pe\_pp}.

\begin{verbatim}
ATOMBASIS
Formamide (aug-cc-pVTZ) and PPs on a water molecule

AtomTypes=6 Charge=0 NoSymmetry Angstrom
Charge=8 Atoms=1 Basis=aug-cc-pVTZ
O     22.931000    21.390000    23.466000
Charge=6 Atoms=1 Basis=aug-cc-pVTZ
C     22.287000    21.712000    22.485000
Charge=7 Atoms=1 Basis=aug-cc-pVTZ
N     22.832000    22.453000    21.486000
Charge=1 Atoms=3 Basis=aug-cc-pVTZ
H     21.242000    21.408000    22.312000
H     23.729000    22.867000    21.735000
H     22.234000    23.026000    20.883000
Charge=8 Atoms=1 Basis=pointcharge ECP=pe_pp
O     21.415000    20.203000    28.463000
Charge=1 Atoms=2 Basis=pointcharge ECP=pe_pp
H     20.840000    20.665000    29.073000
H     22.267000    20.190000    28.899000
\end{verbatim}

\subsection{Polarizable density embedding}\label{subsec:pde}
The polarizable density embedding (PDE) model is an advanced extension to the PE model
that goes beyond purely classical embedding~\cite{pde1,pde2,pde3}.
Instead of resorting to fragment-based multipole expansions of the charge distributions
in the environment, full fragment-based electronic densities are retained, and the
electrostatic interaction is calculated exactly.
Further, an account of non-electrostatic repulsion (also known as Pauli and
exchange repulsion) is included, which is crucial in avoiding electron
spill-out problems.
The PDE model shares the same efficient description of environment polarization as
the PE model, i.e., it uses distributed polarizabilities on environment sites.
It can be used together with the same methods and properties as the PE model,
except molecular gradients and LAOs, which are not available for any method.

A PDE calculation involves two main stages, first a preparation stage, 
which is used to derive the embedding operators and parameters (electrostatic
and repulsion operators; distributed polarizabilities), and, second, the final
embedding calculation.
The preparation stage consists of two parts. Monomer calculations are used 
to obtain the fragment densities of the environment, and these are followed
by dimer calculations to evaluate operator contributions to the core region
from each fragment.
The calculations needed for the setup are somewhat involved, but they are
automated in the \verb|PyFraME| Python package~\cite{pyframe}, and this
is the recommended way of preparing such calculations.
A minimal example input for deriving a PDE potential is given below.
\begin{verbatim}
#!/usr/bin/env python

import pyframe

system = pyframe.MolecularSystem('example.pdb')
core = system.get_fragments_by_number(1)
environment = system.get_remaining_fragments()

system.set_core_region(fragments=core, basis='aug-pcseg-2')
system.add_region(name="environment", fragments=environment,
                    use_multipoles=False,
                    use_polarizabilities=True,
                    use_fragment_densities=True,
                    use_exchange_repulsion=True,
                    basis='loprop-6-31+G*',
                    method='DFT',
                    xcfun='PBE0')

project = pyframe.Project(mpi_procs_per_job=24, jobs_per_node=1, memory_per_job=24000)
project.create_embedding_potential(system)
project.write_potential(system)
project.write_core(system)
\end{verbatim}
Following the setup, the final embedding calculation is run by providing a
\verb|hdf5| file, which is provided by \verb|PyFraME|, with the PDE potential
and activating the PDE option in the \dalinp\ file.
As an example, the following input calculates five excitation energies and
oscillator strengths for a molecule embedded in the potential stored in \verb|final.h5|.
\begin{verbatim}
**DALTON INPUT
.RUN RESPONSE
.DIRECT
.PELIB
*PELIB
.EEF
.PDE
final.h5
**WAVE FUN
.DFT
PBE0
**RESPONSE
*LINEAR
.SINGLE
.DIPLEN
.ROOTS
5
**END OF
\end{verbatim}

\subsubsection{For developers}
Behind the scenes, the setup for a PDE calculations comprises a number of steps,
which are described in the following. Running these manually is tedious and
error-prone, so we recommend using the \verb|PyFraME| Python package for the
setup in general.
For each fragment in the environment, the density needs to be obtained.
These calculations should be run with just the monomer (in vacuum),
and are activated with the \verb|.SAVE DENSITY| option.
\begin{verbatim}
**DALTON
.RUN WAVE FUNCTION
.DIRECT
*PELIB
.SAVE DENSITY
save_density.h5
**WAVE FUNCTIONS
.DFT
PBE0
**END of
\end{verbatim}
In addition to the standard input files, a \verb|hdf5| file
with information about the core region and environment fragment
must be provided. Also, a potential file with the positions of
the polarizable sites must be supplied.
The following fields must be provided in the \verb|hdf5| file:
\begin{verbatim}
HDF5 "save_density.h5" {
GROUP "/" {
   GROUP "core_fragment" {
      DATASET "charges" {
         DATATYPE  H5T_IEEE_F64LE
         DATASPACE  SIMPLE { ( 3 ) / ( 3 ) }
      }
      DATASET "coordinates" {
         DATATYPE  H5T_IEEE_F64LE
         DATASPACE  SIMPLE { ( 3, 3 ) / ( 3, 3 ) }
      }
      DATASET "num_nuclei" {
         DATATYPE  H5T_STD_I64LE
         DATASPACE  SCALAR
      }
   }
   GROUP "fragment" {
      DATASET "charges" {
         DATATYPE  H5T_IEEE_F64LE
         DATASPACE  SIMPLE { ( 3 ) / ( 3 ) }
      }
      DATASET "coordinates" {
         DATATYPE  H5T_IEEE_F64LE
         DATASPACE  SIMPLE { ( 3, 3 ) / ( 3, 3 ) }
      }
      DATASET "num_nuclei" {
         DATATYPE  H5T_STD_I64LE
         DATASPACE  SCALAR
      }
   }
}
}
\end{verbatim}
Note that the dimensions shown here, and in the listings below, depend on the
dimensions (number of atoms / number of basis functions) of the environment fragments
and core region. The \verb|.SAVE DENSITY| calculation will result in a \verb|hdf5| file
to which the following information has been added:
\begin{verbatim}
HDF5 "save_density.h5" {
GROUP "/" {
   GROUP "core_fragment" {
      DATASET "nuclear-electron energy" {
         DATATYPE  H5T_IEEE_F64LE
         DATASPACE  SIMPLE { ( 1 ) / ( 1 ) }
      }
   }
   GROUP "fragment" {
      DATASET "density matrix" {
         DATATYPE  H5T_IEEE_F64LE
         DATASPACE  SIMPLE { ( 28 ) / ( 28 ) }
      }
      DATASET "electric fields" {
         DATATYPE  H5T_IEEE_F64LE
         DATASPACE  SIMPLE { ( 9 ) / ( 9 ) }
      }
      DATASET "energy-weighted density matrix" {
         DATATYPE  H5T_IEEE_F64LE
         DATASPACE  SIMPLE { ( 7, 7 ) / ( 7, 7 ) }
      }
      DATASET "num_bas" {
         DATATYPE  H5T_STD_I32LE
         DATASPACE  SIMPLE { ( 1 ) / ( 1 ) }
      }
      DATASET "num_pols" {
         DATATYPE  H5T_STD_I32LE
         DATASPACE  SIMPLE { ( 1 ) / ( 1 ) }
      }
   }
}
}
\end{verbatim}
The dimer calculations (core region + environment fragment) are activated with the \verb|.TWOINT| option.
The basis functions of the core region should be listed first in the \verb|.mol| file,
followed by the basis functions of the environent fragment.
The \verb|hdf5| output from the preceding monomer calculation should be provided.
\begin{verbatim}
**DALTON
.RUN WAVEFUN
.DIRECT
*PELIB
.TWOINT
save_density.h5
**WAVE FUNCTIONS
.HF
**END OF
\end{verbatim}
This will result add entries for the electrostatic (fragment-density) operator
and exchange-repulsion operator to the output \verb|hdf5| file.
\begin{verbatim}
HDF5 "twoint.h5" {
GROUP "/" {
   GROUP "core_fragment" {
      DATASET "electrostatic matrix" {
         DATATYPE  H5T_IEEE_F64LE
         DATASPACE  SIMPLE { ( 28 ) / ( 28 ) }
      }
      DATASET "exchange-repulsion matrix" {
         DATATYPE  H5T_IEEE_F64LE
         DATASPACE  SIMPLE { ( 28 ) / ( 28 ) }
      }
      DATASET "num_bas" {
         DATATYPE  H5T_STD_I32LE
         DATASPACE  SIMPLE { ( 1 ) / ( 1 ) }
      }
   }
}
\end{verbatim}
From individual dimer calculations on each of the environment fragments,
a final \verb|hdf5| file can be assembled. This file should have the
following structure:
\begin{verbatim}
HDF5 "final.h5" {
GROUP "/" {
   DATASET "electric fields" {
      DATATYPE  H5T_IEEE_F64LE
      DATASPACE  SIMPLE { ( 9 ) / ( 9 ) }
   }
   DATASET "electrostatic matrix" {
      DATATYPE  H5T_IEEE_F64LE
      DATASPACE  SIMPLE { ( 28 ) / ( 28 ) }
   }
   DATASET "exchange-repulsion matrix" {
      DATATYPE  H5T_IEEE_F64LE
      DATASPACE  SIMPLE { ( 28 ) / ( 28 ) }
   }
   DATASET "nuclear charges" {
      DATATYPE  H5T_IEEE_F64LE
      DATASPACE  SIMPLE { ( 3 ) / ( 3 ) }
   }
   DATASET "nuclear coordinates" {
      DATATYPE  H5T_IEEE_F64LE
      DATASPACE  SIMPLE { ( 3, 3 ) / ( 3, 3 ) }
   }
   DATASET "nuclear-electron energy" {
      DATATYPE  H5T_IEEE_F64LE
      DATASPACE  SIMPLE { ( 1 ) / ( 1 ) }
   }
   DATASET "num_bas" {
      DATATYPE  H5T_STD_I32LE
      DATASPACE  SCALAR
   }
   DATASET "num_fields" {
      DATATYPE  H5T_STD_I32LE
      DATASPACE  SCALAR
   }
   DATASET "num_nuclei" {
      DATATYPE  H5T_STD_I32LE
      DATASPACE  SCALAR
   }
}
}
\end{verbatim}

\subsubsection{London orbitals with PDE}
The PDE model supports the calculation of NMR shieldings using London orbitals.
Additional magnetic gradient terms need to be evaulated when constructing the
embedding potential. The keyword \verb|.TWOLAO| should be activated under the
\verb|*PELIB| section (see input listing below), which will request the
generation of the \verb|london electrostatic matrix| and
\verb|london exchange-repulsion matrix| in the \verb|hdf5| file.

\begin{verbatim}
**DALTON
.RUN WAVEFUN
.DIRECT
*PELIB
.TWOINT
save_density.h5
**WAVE FUNCTIONS
.HF
**END OF
\end{verbatim}


\subsection{Fast Multipole Method}\label{subsec:FMM}
The computational scaling of the PE model is quadratic with the number of
polarizable sites. For large systems, this can become the bottleneck of PE
calculations. To overcome this, fast summation schemes, such as the fast
multipole method (FMM) can be used (see Ref. \citenum{scheurer2021efficient}).
The FMM implementation in PElib is used to evaluate both the static fields
from multipoles as well as the fields from induced dipoles. In practice, FMM
is beneficial (faster than traditional direct summation) for systems exceeding
about 8000 environment sites. FMM can be activated with:

\begin{verbatim}
**DALTON
.RUN WAVEFUN
.DIRECT
*PELIB
.FMM
**WAVE FUNCTIONS
.HF
**END OF
\end{verbatim}

The accuracy of FMM is controlled by two parameters, the multipole acceptance criterion, $\theta$,
and the multipole expansion order. The multipole acceptance criterion determines whether two
boxes, $A$ and $B$ can be considered far separated, according to
\begin{equation}
R\theta > R_A + R_B.
\end{equation}
Smaller values of $\theta$ are more accurate since more of the interactions are placed in the
short range, but are also more expensive. The limit of $\theta=0$ corresponds formally to a
direct summation. Typical values of $\theta$ are around $0.5$. The expansion order controls
to which order the multipole expansion in FMM is taken. Higher orders are more accurate, but
also more expensive. The FMM implementation is based on an adaptive octree code, which
boxifies the system such that boxes contain at most NCRIT particles. Smaller values of NCRIT
will result in more boxes. The acceptance criterion ($\theta$), expansion order, and critical
number of particles can be changed as:

\begin{verbatim}
**DALTON
.RUN WAVEFUN
.DIRECT
*PELIB
.FMM
.THETA
0.40       ! default
.EXPANSION ORDER
6          ! default
.NCRIT
64         ! default
**WAVE FUNCTIONS
.HF
**END OF
\end{verbatim}

\subsection{Approximate environment coupling}
After applying FMM to treat the interaction between induced dipoles, all
classical and quantum-classical interaction parts of a PE calculation scale
linearly with respect to the number of environment sites. For some systems,
the evaluation of electric fields from the electronic density in the quantum
part and the corresponding evaluation of Fock matrix elements due to the
induced dipoles can become a dominating factor in the cost of a PE calculation.
Various approximations to the interaction between the quantum region and the
environment can be introduced to reduce the cost, with only a very minor impact
on the quality of the calculation (see Ref. \citenum{reinholdt2021fast}). One
option is to use a single-center multipole expansion. Instead of using the full
interaction between the quantum region and the environment, electric fields
beyond some cutoff radius are calculated via a single-center multipole expansion.
The expansion is generally only convergent beyond some radius, and relatively
high expansion orders should be used for good accuracy. Below is an example of
a single-center expansion with a cutoff radius of 40 bohr and an interaction
multipole expansion order of 8.

\begin{verbatim}
**DALTON INPUT
.RUN WAVEFUN
.DIRECT
.PEQM
*PEQM
.FMM
.THETA
0.4
.EXPANSION_ORDER
6
.INT_SCHEME
SINGLE_CENTER
.INT_ORDER
8
.INT_RCUT
40.0
**WAVE FUNCTIONS
.DFT
CAMB3LYP
**END OF
\end{verbatim}

Another available option is to use a multi-center multipole expansion with
ESP-fitted multipoles. The multi-center expansion improves the short-range
convergence of the expansion. Expansion orders should generally be kept low.
Below is an example of an ESPF expansion with up to quadrupoles and a cutoff
radius of 20 bohr. Note the difference between `INT\_RCUT` (determines where
the approximate coupling can be applied) and `INT\_RFIT` (determines the
radius within which solvent sites are used as grid points for the ESP fitting).

\begin{verbatim}
**DALTON INPUT
.RUN WAVEFUN
.DIRECT
.PEQM
*PEQM
.FMM
.THETA
0.4
.EXPANSION_ORDER
6
.INT_SCHEME
ESPF
.INT_ORDER
2
.INT_RCUT
20.0
.ESPF_GRID
SOLVENT
.INT_RFIT
20.0
**WAVE FUNCTIONS
.DFT
CAMB3LYP
**END OF
\end{verbatim}

Another option is to create local expansions of the electronic electric fields
in the environment, similar to FMM. Here, the near-field/far-field cutoff is
controlled by a multipole acceptance criterion, $\theta$. Smaller values of
$\theta$ (`INT\_THETA`) place more interactions in the near-field, and are thus
more accurate but also more expensive. Below is an example of this approach,
with $\theta=0.3$ and a local expansion up to 7th order (`INT\_ORDER`).

\begin{verbatim}
**DALTON INPUT
.RUN WAVEFUN
.DIRECT
.PEQM
*PEQM
.FMM
.THETA
0.4
.EXPANSION_ORDER
6
.INT_SCHEME
FMM_FOCK
.INT_ORDER
7
.INT_THETA
0.3
**WAVE FUNCTIONS
.DFT
CAMB3LYP
**END OF
\end{verbatim}
